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Abstract 

The interaction of polymers with small-scale velocity gradients can trigger a coil-stretch transition in the polymers. 
We analyze this transition within a direct numerical simulation of shear turbulence with an Oldroyd-B model 
for the polymer. In the coiled state the lengths of polymers are distributed algebraically with an exponent a — 
2y — 1/De, where 7 is a characteristic stretching rate of the flow and De the Deborah number. In the stretched 
state we demonstrate that the length distribution of the polymers is limited by the feedback to the flow. 
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1. Introduction 

Dilute solutions of long polymers show a fasci- 
nating variety of unusal flow phenomena [1,2]. Re- 
cent effects include the formation of vortex pairs in 
viscoelastic Taylor- Couette flow [3-5] and a phe- 
nomenon called elastic turbulence [6] . While many 
phenomena in the laminar regime can be under- 
stood within simple constitutive equations [1,2], 
the interplay with turbulence has remained less 
clear. In particular, the drag reduction by minute 
amounts of long polymers in turbulent flows [7] 
awaits explanation. Various experiments have es- 
tablished that drag reduction comes about when 
the relaxation time of the polymers is comparable 
to the fluctuation times in the velocity field [8,9]. 
And it is known that a long polymer undergoes 
a transition to an uncoiled state when exposed to 
sufficient strain [10,11]. Our aim here is to study 
the interaction of a polymer with a turbulently 



fluctuating flow, in particular the transition from 
the coiled to the stretched state as the internal re- 
laxation rate becomes comparable to the velocity 
gradients, and to compare to the theoretical pre- 
dictions of Balkovsky et al. [12]. Specifically, they 
predict an algebraic distribution for the trace of 
the configuration tensor with an exponent that de- 
pends linearly on the Deborah number and an in- 
ternal stretching rate of the fiow. Moreover, if the 
strain is so strong as to uncoil the fiow, this process 
comes to a stop through its effect on small-scale 
fiuctuations of the flow. 



2. Model and numerical implementation 

We model the polymer by an upper convected 
Maxwell fluid with a single relaxation time A [13]. 
The polymer is characterized by the conformation 
tensor Cy . Within a dumb-bell model [1] this tensor 
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can be formed from the end-to-end distance vector 
R as Cij — {RiRj) where (■ • ■ ) denotes a thermal 
average. In the coiled state the polymers are spher- 
ical and c can be normalized so that c^j* = 5ij and 
Trc = 3. In the Maxwell model the relaxation of 
the polymers is assumed to be linear, so that the 
dynamical equations become 



D 
Dt 



Cij - CikdkVj - CjkdkVi = 



X 



(1) 



with A as the time constant and D/ Dt as the con- 
vective derivative. The polymers are assumed to 
be diluted in a Newtonian solvent so that for the 
flow properties also the Newtonian shear viscosity 
has to be included. This then defines an Oldroyd- 
B model [14]. The Oldroyd-B model is able to de- 
scribe two main features of viscoelasticity, namely 
normal stress differences and stress relaxation, and 
has been succesfuUy applied in the study of vortex 
pairs in Couette- Taylor flow [5] . The main caveat 
of this model, the divergence of the conformation 
tensor when the shear rate of the flow exceeds 1/ A 
is usually overcome with the FENE model [1,15]. 
However, as predicted by Balkovsky ct al. [12] and 
demonstrated below the extension of the polymer 
is stabilized through the feedback on the flow once 
it begins to expand. 

We want to study the turbulence of a flow 
bounded by two surfaces and driven by a vol- 
ume force that maintains a constant mean shear 
gradient S (see Fig. 1). This gives an external 
length scale L (distance between surfaces) and 
a velocity scale U (twice the velocity difference 
across the gap).[£| The Reynolds number then is 
Re = UL/v with v being the kinematic viscosity, 
and the Deborah number (characterizing the poly- 
mer relaxation) is De — XU/L. The dimensionless 
Navier-Stokes and Oldroyd-B equations become 

1 



dtVi + {vkdk)vi = -dtP + -^djjVi + d^rf^ + f, 

On - 5, 



(2) 



De 



(3) 



^ The unusual appearance of a factor 2 is connected with 
the dimensionless shear rate S = 0.5 used in the code. 
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Fig. 1. Model geometry: The flow is restricted to the volume 
between the shaded planes, with no-flux, free-slip boundary 
conditions at the planes, ft is periodically continued in the 
y- and 2-direction. 

The relation of stress and conformation tensor 
is given by 



DeRe 



(cij Sij) . 



(4) 



The parameter s in (4) models the strength of the 
feedback of the polymers on the turbulent flow. 
Together with the condition of incompressibility, 



dtV, = , 



(5) 



and a prescribed external force fi , the set of equa- 
tions (2) - (4) form a closed system for the vari- 
ables Vi, Cij and p. 

The advantage of writing the Oldroyd-B equa- 
tion in terms of the conformation tensor c instead 
of the polymer stress tensor is that even for s — 
0, when the polymer has no influence on the veloc- 
ity field and rf^ — 0, it is still possible to monitor 
c and thus the stretching and rotation of the poly- 
mer by the flow. 

The model geometry (see Fig. 1) resembles pla- 
nar Couette flow, i.e. flow between infinite, paral- 
lel plates, but with free-slip (or stress- free) bound- 
ary conditions along the planes instead of no-slip 
boundary conditions. With these boundary con- 
ditions all fields can be expanded in Fourier se- 
ries and efhcient algorithms can be used. In our 
dimensionless units the box has dimensions Q — 
[0,1] X [0,Ly] X [0,L,]. 

While any quantity is taken to be periodic in 
the y- and z-direction, boundary conditions in the 
x-direction are somewhat more complicated and 
must be formulated for each quantity separately. 
For the velocity, they are 
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dxVy\x=0.,l = 

dxVz\x=a.i - 



--0; 

0, 
0. 



(no flux) 
(no stress) 



(6) 



This translates for the solvent stress tensor rf- into 



Tij \x=o,i = for ij = xy, xz, yx, zx, 
dxT^j \x=o,i = for ij = XX, yy, zz, yz, zy. 



(7) 



The boundary conditions on r?- and Ctj must be 
the same as on r?- , 



^ij la:— 0,1 



for ij = xy, xz, yx, zx. 



dxCtj\x=o.i = for ij 



(8) 

xx,yy,zz,yz,zy. 

Instead of considering the quantities Vi and Cy 
over fl and with boundary conditions (6) and (8) , 
one can extend them to 2f7 := [—1, 1] x [0, Lj,] x 
[0, Ly] by reflection at the x = 0-plane. The equa- 
tions of motion (2) and (3) are invariant under 
this transformation. Along the boimdary of 2Q, 
any quantity then obeys periodic boundary condi- 
tions and can be expanded in a three-dimensional 
Fourier scries. Mirror symmetry and boundary 
conditions restrict modes in the a;-direction to 
either sines or cosines. 

The simulation code is actually written for 
mixed Fourier {y, z) and sine/cosine {x) expan- 
sion coefficients. The reahty condition, relations 
between coefficients resulting from incompress- 
ibility and de-aliasing by the 2/3 rule [16] lead to 
a further reduction in the number of independent 
variables. In the end one is left with 

{4. + mx) + {^ + mx){N, + Ny + NyN,) (9) 

independent coefficients, where Ni = 2iVj/3 for a 
spatial grid of (A^ + 1) x Ny x A^^ points. 

The numerical algorithm employed for time inte- 
gration is a 5(4)-Dormand-Prince scheme [17] with 
flxed step size. This scheme uses five evaluations of 
the right hand side to calculate an approximation 
foT:x{t + At), 
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x{t + At) = J2 + 0{Af) , (10) 

i=0 

where Kj = F{xj,tj) , 

i-i 

Xj = x{t) -I- Ai ^ ajiKi , 



j=0 



with coefficients aij and 7, chosen to minimize the 
error accumulated over At [17]. Then, using the 
AevwaXwe F{x{t + At),t + At) « ^(xe.ie) (which 
can be kept for the next integration step), a second 
approximation 

6 

x'{t + At) = Y,<iKi + 0{At^) (11) 

of lower order can be obtained. The difference 5i = 
Xi — x'i is used as an estimate for the integration 
error. In standard applications this error is then 
used to adapt the step size. However, in a turbulent 
ffow with its ever present fluctuations, variations 
in this step size are small. So instead we monitor 
this error estimate by calculating a combination of 
relative and absolute error, 



Err 



max ■ 



\6i 



i |a;i|-|-e 
10~^ an accuracy of Err < 10~^ 



(12) 



was 



With e 

maintained throughout the integration. 

The volume force fi was chosen so as to maintain 
the mean shear profile (dimensionless form) 



V'yix) 
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M 

E 

m=l 



COS {(2ni — l')n.i:) 
(2m- 1)2 



(13) 



tj = t + jjAt, 



as an approximation to a linear profile Vy{x) = 
S{x — 1/2). In the simulations, M = 5, and the 
Fourier components contained in this mean profile 
were kept fixed. Thus the force is dynamically ad- 
justed to maintain the prescribed mean shear pro- 
file. 

One way to initialize the program is to use a ran- 
dom velocity perturbation on the laminar profile 
and the conformation tensor for the laminar pro- 
file. In order to obtain a turbulent state the pertur- 
bation has to be about as large as the laminar pro- 
file. In turns out, however, that even for moderate 
Deborah numbers De w 0.1, . . . ,1 this leads to an 
extreme overshooting of the conformation tensor 
and a breakdown of the integration routine. 

As a remedy to this we used initial conditions 
from a turbulent Newtonian fiow with equilibrium 
configuration tensor, = 6ij. The system was 
then allowed to relax to a statistically stationary 
turbulent state for a small value of the Deborah 
number. One of these states was then taken as ini- 
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Fig. 2. Probability density functions (PDFs) of Trc in a 
turbulent flow at different De. There are hardly any regions 
with negative Tr c for De = 0.5 (open triangles) and a sig- 
nificant fraction for De = 0.7 (open circles). Adding an ar- 
tifical diffusivity of k = 1 /2000 removes the negative values 
and leaves the positive distribution essentially unchanged 
(filled circles). Parameters: Re = 1000, S = dx{vy) = 0.5 
and s = 0. Geometry: Ly = in, = 27r, = 32, 
Ny = 128, = 64. The PDFs are averaged over 10 snap- 
shots taken at time intervals At = 5. 

tial condition for a simulation at a higher De. Pro- 
ceeding in this manner at each change of De the 
system relaxed to a statistically stationary state 
within a few time steps and no further numerical 
integration problems were encountered. 

Another numerical problem for Oldroyd-B sim- 
ulations is connected with the fact that the con- 
formation tensor Cij is by definition positive def- 
inite, and that this property is preserved by the 
Oldroyd-B equations [2] but not in the numerical 
representation. Through the various Fourier trans- 
forms and truncations, this positivity can be lost. 
If that happens the Oldroyd-B model can develop 
Hadamard instabilities, and lose its evolutionary 
character [18]. Thus, if regions with negative c ap- 
pear, numerical instabilities can arise and the in- 
tegration can blow up [19]. Fig. 2 shows that such 
regions indeed do accur, but remain limited to a 
small fraction of the whole domain, provided De is 
not too large. 

As a partial cure to this problem, Beris et al. [19- 
21] introduce an artificial stress diffusion — kAc in 
the equation of motion for the conformation ten- 
sor, where n can bc' chosen so that it stabilizes the 
numerics but has negligible influence on the re- 
sults. Although for the small values of De in Fig. 2 
no numerical instability seems to occur, the intro- 
duction of a small k does reduce the probability 
for negative Tr c, while the impact on positive val- 



ues is much less significant. In the following, cal- 
culations for De < 0.5 are without artificial stress 
diffusivity (k = 0). The number of points with 
Tr c < is then negligible. For De > 0.5, a small 
K = 1 /2000 is employed to reduce negative tails of 
PDFs, and for Z? = 2.0 it is actually necessary to 
set K = 1/1000 to avoid numerical instability. In 
both cases, a certain probability for negative Trc 
is still unavoidable. 

The simulation code is written in Fortran 90 (ex- 
cept the special FFT routines for Cray T90 hard- 
ware) . It is based on the code used previously [22] 
for the same model without polymer. Most simu- 
lations were done on a Cray T90 vector machine at 
the John von Neumann Institute for Computing at 
the Research Centre Jiilich. A typical rmi with a 
resolution of TV^,; = 32, TV^ = 128, A^^ = 64 over 50 
time units costs roughly 26500 CPU seconds and 
needs 62 megawords of memory. 



3. Statistics of polymer elongation 

Balkovsky et al. [12] discuss the behaviour of the 
PDF for the polymer extension R = (Tr c)^/^. For 
a passive polymer, they predict a power law decay 
of the tails {R oo), as long as the (dimension- 
less) relaxation time of the polymer De is smaller 
than the inverse shear rate 1/(27) of the underly- 
ing turbulent flow, 

V{R) - i?r"-i . (14) 

The exponent a depends on De and 7. For relax- 
ation times near the inverse Lyapunov exponent, 
a is given in linear approximation by 



For Deborah numbers smaller than Deer = 1/ (27) 
the exponent a is positive and (14) actually de- 
scribes an algebraic decay. This corresponds to the 
majority of the polymer molecules being near their 
equilibrium elongation. However, for values above 
the critical Deborah number Deer = 1/(27) the 
exponent a becomes negative and the probability 
density function for R is no longer normalizable. 
At this point the assumption of a passive polymer 
becomes unphysical. Balkovsky et al. [12] then ar- 
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Fig. 3. PDFs of Tr c in turbulent flow for different values of 
De. The straight lines in the lower frame follow from alge- 
braic fits over the ranges over which they are shown. Every 
data point corresponds to one bin. Parameters; Re = 1000, 
da:{vy) = 0.5 and s = 0. k = 1/2000 for De > 0.6. Geome- 
try: Ly = 4tt, = 2tt, Nx = 32, Ny = 128, = 64. The 
initial condition was the same for all De. The PDFs are 
averaged over 10 snapshots taken at time intervals At = 5. 

gue that the feedback from the polymer to the flow 
field limits the extension, such that V{R) decays 
again for large R. The majority of the molecules is 
then stretched to some finite value of R, at which 
V{R) is maximal. The contribution of the poly- 
mer stress to the right hand side of the Navier- 
Stokes equation should then be of the same order 
of magnitude as the advective term v ■ Vv. 

The quantity directly accessible from the simu- 
lations is Trc = The PDFs for i?^ and R are 
related through 

V{R)dR = V{R^)dR^ = 2RV{R^)dR. (16) 

A power law distribution for one quantity implies 
one for the other, and the exponents are related by 

P(Tr c) ~ (Tr c)"'^ <^ V{R) ~ R'"'^ , (17) 

with a = 2{0 - 1). 

Fig. 3 demonstrates these power tails in PDFs 
obtained from turbulent simulations for different 




Fig. 4. Exponents extracted from Fig. 3 {De < 0.8, squares) 
and the equivalent for (De = 2.0, s = 0.01, open circle). 
The regression line derived from the solid squares intersects 
with a = for 1/ De = 0.61 = 27. Note that the open circle 
for the situation of uncoiled polymers also lies on the line. 

De. The exponents (3 are obtained from the slope 
in a log- log plot (lower graph). The region where 
the power law actually holds starts when 'P(Trc) 
drops below about 10~^ and extends right up to 
the statistical accuracy limit for larger Tr c. 

From the exponents determined in Fig. 3, the 
Lyapunov exponent 7 (and with it, the critical 
Deborah number) can be obtained by plotting a 
as a function oil/De, as in Fig. 4, and extrapolat- 
ing the linear relationship to a = 0. Note that the 
three points corresponding to De = 0.6, 0.7 and 
0.8 fit the line surprisingly well. The deviation for 
De = 0.5 could be attributed to the fact that only 
this simulation is done with k = 0. Also, De = 0.5 
is farthest from the (extrapolated) critical Debo- 
rah number. 

From the data on the exponents one can extract 
that for De larger than 1/0.61 the polymer dis- 
tribution without feedback is no longer normaliz- 
able, the polymers go into a stretched state. This 
stretching is terminated by their feedback on the 
flow, as we will now demonstrate. To this end we 
analyze the time evolution of (Tr c) in Fig. 5. The 
trace starts at the equilibrium value Tr c = 3 at t = 
0. For De < De^r it increases and then stabilizes 
at a finite value, indicating a certain stretching of 
the polymer. For De > Deer it increases monoton- 
ically. if there is no feedback to the flow (s = 0), 
but stabilizes at finite values for s ^ 0. The larger 
the feedback, the smaller the mean values. 
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Fig. 5. Time evolution of (Trc) for De > Dccv With- 
out feedback through the flow the polymer length diverges. 

With feedback (s ^ 0) it stabilizes. The horizontal dot- 
ted lines mark the time average of the solid points; these 
averages decrease with increasing coupling strength s. 

In an additional calculation for De = 2.0 and 
s = (see Fig. 5), (Trc) grows more or less expo- 
nentially, as expected, with a rate dln(Tr c)/(it « 
0.18. The difference between l/Dccr = 0.61 and 
1/De = 0.5 gives a growth rate of about 0.1. The 
difference is perhaps connected to the use of differ- 
ent averages in the calculation of the exponents. 



4. Conclusions 

The full numerical simulations for an Oldroyd- 
B constitutive equation in a turbulent shear flow 
support the calculations of Balkovsky et al. and 
demonstrate and algebraic distribution of polymer 
lengths below the coil-stretch transition. Above the 
transition the calculations also show that the po- 
tentially infinite stretching of the polymers is lim- 
ited by their modulation of the small-scale flow 
properties. 
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